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ABSTRACT 


We study the effects of planetary late migration on the gas giants obliqui¬ 
ties. We consider the planetary instability models from [Nesvorny fc Morbidelli 


(2012), in which the obliquities of Jupiter and Saturn can be excited when the 
spin-orbit resonances occur. The most notable resonances occur when the Sj and 
sg frequencies, changing as a result of planetary migration, become commensu¬ 
rate with the precession frequencies of Jupiter’s and Saturn’s spin vectors. We 
show that Jupiter may have obtained its present obliquity by crossing of the sg 
resonance. This would set strict constrains on the character of migration during 
the early stage. Additional effects on Jupiter’s obliquity are expected during the 
last gasp of migration when the sy resonance was approached. The magnitude 
of these effects depends on the precise value of the Jupiter’s precession constant. 
Saturn’s large obliquity was likely excited by capture into the sg resonance. This 
probably happened during the late stage of planetary migration when the evo¬ 
lution of the Sg frequency was very slow, and the conditions for capture into the 
spin-orbit resonance with sg were satisfied. However, whether or not Saturn is 
in the spin-orbit resonance with sg at the present time is not clear, because the 
existing observations of Saturn’s spin precession and internal structure models 
have significant uncertainties. 
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Subject headings: celestial mechanics — planets and satellites: individual (Jupiter, 
Saturn) — planets and satellites: dynamical evolution and stability 


1. Introduction 


It is believed that the orbital architecture of the solar system was significantly altered 
from its initial state after the dissipation of the protosolar nebula. The present architecture 
is probably a result of complex dynamical interaction between planets, and between planets 
and planetesimals left behind by planet formation. This becomes apparent because much of 
what we see in the solar system today can be explained if planets radially migrated, and/or if 


they evolved through a dynamical instability and reconfigured to a new state (e.g., Malhotra 


1995 

Thommes et ah 

1999 

Tsiganis et ah 

2005 


While details of this process are not known exactly, much has been learned about it over 
decade by testing various migration/instability models against various constraints. Some of 
the most important constraints are provided by the terrestrial planets and the populations of 


small bodies in the asteroid and Kuiper belts (e.g.. 

Gomes et ah 

2005 

Minton & Malhotra 

2009 

2011 

Levison et ah 2008 Morbidelli et ah 

2010; 

Nesvorny 

2015 

). Processes 


related to the giant planet instability/migration were also used to explain capture and orbital 


distribution of Jupiter and Neptune Trojans (e.g.. 

Morbidelli et ah 2005 

Nesvorny & 

Vokrouhlicky 

2009; 

Nesvorny et ah 

2014a|) and irregular satellites (e.g.. 

Nesvorny et ah 


2014b). Some of the most successful instability/migration models developed so far 


postulate that the outer solar system contained additional ice giant that was ejected into 


interstellar space by Jupiter (e.g., Nesvorny 2011; Nesvorny & Morbidelli 2012, Batygin 


et ah 2012). The orbits of the four surviving giant planets evolved in this model by 
planetesimal-driven migration and by scattering encounters with the ejected planet. In this 
work, we use this framework to investigate the behavior of Jupiter’s and Saturn’s obliquities. 


The obliquity, e, is the angle between the spin axis of a planet and the normal to its 
orbital plane. The core accretion theory applied to Jupiter and Saturn implies that their 
primordial obliquities should be very small. This is because the angular momentum of the 
rotation of these planets is contained almost entirely in their massive hydrogen and helium 
envelopes. The stochastic accretion of solid cores should therefore be irrelevant for their 
current obliquity values (see Lissauer & Safronov 1991 for a discussion), and a symmetric 
inflow of gas on forming planets should lead to £ ~ 0. The present obliquity of Jupiter is 
ej = 3.1°, which is small, but not quite small enough for these expectations, but that of 
Saturn is £s = 26.7°, which is clearly not. 
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Ward & Hamilton (2004) and Hamilton & Ward (2004) noted that the precession 
frequency of Saturn’s spin axis has a value close to Sg — —0.692 arcsec yr“^, where Sg is the 
mean nodal regression of Neptune’s orbit (or, equivalently, the eighth nodal eigenfrequency 


of the planetary system; e.g., Applegate et ah 1986; Laskar 1988). Similarly, Ward & 


Canup (2006) pointed out that the precession frequency of Jupiter’s spin axis has a value 
close to Sj ~ —2.985 arcsec yr“^, where Sj is the mean nodal regression of Uranus’s orbit. 
These findings are significant because they raise the possibility that the current obliquities 
of Jupiter and Saturn have something to do with the precession of the giant planet orbits. 


Specihcally, Ward & Hamilton (2004) and Hamilton & Ward (2004) suggested that the 


present value of Saturn’s obliquity can be explained by capture of Saturn’s spin vector in 
a resonance with sg. They proposed that the capture occurred when Saturn’s spin vector 
precession increased as a result of Saturn’s cooling and contraction, or because Sg decreased 
during the depletion of the primordial Kuiper belt. They showed that, if the post-capture 
evolution is conveniently slow, the spin-orbit resonance (also known as the Cassini resonance, 
see Section 2) is capable of increasing Saturn’s obliquity to its current value. 

While changes of precession or sg during the earliest epochs could have been important, 
it seems more likely that capture in the spin-orbit resonance occurred later, probably as a 
result of planetary migration (Boue et ah 2009). This is because both sr and sg signih- 


cantly change during the instability and subsequent migration. Therefore, if the spin-orbit 
resonances had been established earlier, they would not survive to the present time. |Boue 


et ah (2009) studied various scenarios for resonant tilting of Saturn’s spin axis during 


the planetary migration and found that the present obliquity of Saturn can be explained 
only if the characteristic migration time scale was long and/or if Neptune’s inclination was 
substantially excited during the instability. Since Neptune inclination is never large in the 
instability/migration models of Nesvorny & Morbidelli (2012; hereafter NM12), typically 


zn < 1°, the migration timescales presumably need (note that Boue et ah (2009) did not 


investigated these low-z cases in detail) to be very long (see Section 3.3). Interestingly, these 


very long migration timescales are also required from other constraints (e.g., Morbidelli et 


al. 2014; Nesvorny 2015). They could be achieved in the Nesvorny & Morbidelli (2012) 


models if Neptune interacted with an already depleted planetesimal disk during the very 


last stages of the migration process. As for the obliquity of Jupiter, Ward & Canup (2006) 


suggested that the present value is due to the proximity of the spin precession rate to the sj 
frequency. 


In fact, the obliquities of Jupiter and Saturn represent a much stronger constraint on 
the instability/migration models than was realized before. This is because the constraints 
from the present obliquities of Jupiter and Saturn must be satisfied simultaneously (McNeil 
&: Lee 2010; Brasser & Lee 2015). For example, in the initial conhguration of planets in the 
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NM12 models, the Sg frequency is much faster than the precession constants of both Jupiter 
and Saturn. This means that the Sg mode, before approaching Saturn’s precession frequency 
and exciting its obliquity, must also have evolved over the precession frequency of Jupiter’s 
spin vector. This leads to a conundrum, because if the crossing were slow, Jupiter’s obliquity 
would increase as a result of capture into the spin-orbit resonance with sg. If, on the other 
hand, the general evolution at all stages were fast, the conditions for capture of Saturn into 
the spin-orbit resonance with sg may not be be met (e.g., 
obliquity would stay smallj^ 


Bone et al. 2009), and Saturn’s 


A potential solution of this problem would be to invoke fast evolution of sg early on, 
during the crossing of Jupiter’s precession frequency, and slow evolution of Sg later on, such 
that Saturn’s spin vector can be captured into the spin-orbit resonance with sg during the 
late stages. This can be achieved, for example, if the migration of the outer planets was faster 
before the instability, and slowed down later, as the outer Solar System progressed toward 
a more relaxed state. As we show in Section 3.1, the jumping-Jupiter models developed in 
NM12 provide a natural quantitative framework to study this possibility. 


In Sec. we hrst briefly review the general equations for the spin-orbit dynamics. Then, 
in Sec. we investigate the behavior of the spin vectors of Jupiter and Saturn in the 
instability/migration models of NM12. We hnd that the constraints posed by Jupiter’s 
and Saturn’s obliquities can be satished simultaneously in this class of models, and derive 
detailed conditions on the migration timescales and precession constants that would provide 
a consistent solution. 


2. Methods 

2.1. Parametrization using non-singular spin vector components 

Consider a planet revolving about the Sun and rotating with angular velocity oj about 
an instantaneous spin axis characterized by a unit vector s. With solar gravitational torques 
applied on the planet, u remains constant, but s evolves in the inertial space according to 


^Note that the present obliquities of Uranus and Neptune are not an important constraint on planetary 
migration, because their spin precession rates are much slower than any secular eigenfrequencies of orbits 
(now or in the past). Therefore, the secular spin-orbit resonances should not occur for these planets, and 
giant impacts may be required to explain their obliquities (e.g., Morbidelli et al. 2012 and references 
therein). 
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(e.g., Colombo 1966) 


(jc 

— = —a (c ■ s) (c X s) 
dt V M y 


where c denotes a unit vector along the orbital angular momentum, and 

3 GMq J 2 Q 
2 ujb^ X + I 


( 1 ) 

( 2 ) 


is the precession constant of the planet. Here, G is the gravitational constant, Mq is the 
mass of the Sun, b = ay/l — e^, where a and e are the orbital semimajor axis and ec¬ 
centricity, J 2 is the quadrupole coefficient of the planetary gravity held, and A = G/MR^ 
is the planetary moment of inertia G normalized by a standard factor MR^, where M is 
planet’s mass and R its reference radius (to be also used in the dehnition of J 2 ). The term 
q = \ is an effective, long-term contribution to the quadrupole coefficient 

due to the massive, close-in regular satellites with masses rrij and planetocentric distances 
Oj, and I = '^j^mjOjrij)/{MR^ oj) is the angular momentum content of the satellite orbits 
{rij denotes their planetocentric mean motion) normalized to the characteristic value of the 
planetary rotational angular momentum. 


For both Jupiter and Saturn, q is slightly dominant over J 2 in the numerator of the 
last term in Eq. (|^, while I is negligible in comparison to A in the denominator. Spacecraft 
observations have been used to accurately determine the values of J 2 , q and 1. On the 
other hand, A cannot be derived in a straightforward manner from observations, because it 


depends on the structure of planetary interior. Using models of planetary interior. Helled 


et ah (2011) determined that Jupiter’s Aj is somewhere in the range between 0.2513 and 


0.2529 (here rescaled for the equatorial radius R = 71,492 km of the planet). This would 
imply Jupiter’s precession constant aj to be in the range between 2.754 arcsec yr“^ and 
2.772 arcsec yr“^. These values are smaller than the ones considered in Ward Sz Canup 


(2006), if their proposed small angular distance between Jupiter’s pole and Cassini state C 2 
is correct. 


Similarly, Helled et al. 


(2009) determined Saturn’s precession constant as to be in the 
range between 0.8443 arcsec yr“^ and 0.8447 arcsec yr“^. Again, these values differ from 
those inferred by Ward & Hamilton (2004) and [Hamilton fc Ward (2004), if a resonant 
conhnement of the Saturn’s spin axis in their proposed scenario is true. This difference has 
been noted and discussed in Boue et al. (2009). If Helled’s values are correct, Saturn’s spin 
axis cannot be presently locked in the resonance with sg. However, it seems possible that 


of as derived in Helled et al. (2009) may have somewhat larger uncertainty than reflected 


by the formal range of the inferred A values. Also, the interpretation is complicated by the 
past orbital evolution of planetary satellites which may have also contributed to changes of 
as (and aj). For these reason, and in the spirit of previous studies (e.g.. Ward & Hamilton 
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2004 Ward & Canup 2006; Boue et al. 2009), here we consider a wider range of the 


precession constant valnes a for both Jnpiter and Satnrn. 

Assnming a constant for a moment, a difficnlt element preventing a simple solntion 
of Eq. ([^ is the time evolntion of c. This is becanse mntnal planetary interactions make 
their orbits process in space on a characteristic timescale of tens of thonsands of years and 
longer. In addition, dnring the early phase of planetary evolntion, the precession rates 
of planetary orbits may have been faster dne to the gravitational torqnes from a massive 
popnlation of planetesimals in the trans-Neptnnian disk. In the Keplerian parametrization 
of orbits, the nnit vector c depends on the inclination / and longitnde of node snch that 
= (sin / sin — sin / cos ff, cos /). Traditionally, the difhcnlties with the time dependence 
in Eq. ([^ are resolved nsing a transformation to a reference frame hxed on an orbit, where 
c^ = (0,0,1). 

The transformation from the inertial to orbit coordinate frames is achieved by applying 
a 3-1-3 rotation seqnence with the Enlerian angles (f2,J, —f2). This transforms Eq. ([^ to 
the following form 

(jc 

— = - [a (c ■ s) c + h] X s, (3) 

where now the planetary spin vector s is expressed with respect to the orbit frame, and c is 
now a nnit vector along the z-axis. In effect, the time dependence has been moved to the 
vector qnantity = {A, B, —2C) with 


A = cos QI — sin I sin fl, 
B = sin fl I + sin I cos 
C = sin^ I/2Q, 


(4a) 

(4b) 

(4c) 


and the over-dots mean time derivative. 


A farther development consists in introdncing complex and non-singnlar parameter 
( = sin(//2) exp(zfl) [t = ^/—l is the complex nnit) that replaces I and fl. First order 
pertnrbation theory for qnasi-circnlar and near-coplanar orbits indicates ( for each planet 
can be expressed nsing a hnite nnmber of the Fonrier terms with the s* freqnencies nniqnely 
dependent on the orbital semimajor axes and masses of the planets, and amplitndes set 
by the initial conditions (e.g., Bronwer & Clemence 1961). In the models from non¬ 


linear theories or nnmerical integrations, ( can still be represented by the Fonrier expansion 


( = ^ Ajexp[*(sit -|- 0j)] (e.g., Applegate et al. 1 1986 Laskar 1988), with the linear terms 
having typically the largest amplitndes Aj. 

As discnssed in Section 1, terms with present freqnencies sy ~ —2.985 arcsec yr“^ and 
Sg ~ —0.692 arcsec yr“^ are of a particnlar importance in this work. The sy-term is the largest 
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in Uranus’ ( representation, and the Sg-term is the largest in Neptune’s ( representation, 
and these terms also appear, though with smaller amplitudes, in the ( variable of Jupiter 
and Saturn, because mutual planetary interactions enforce all fundamental frequency terms 
to appear in all planetary orbits. In terms of (, Eqs. (4a-c) become 


A + iB 
C 


\/i-CC 


2i \ at 




(5a) 

(5b) 


where ( is complex conjugate to (. For small inclinations, relevant to this work, we therefore 
hnd that ^ ~ 2 {d(/dt), and C ~ 0. 


Another important aspect is of the problem is that Eq. ([^ derives from a Hamiltonian 

= |(c-s)^ + h-s, (6) 

such that 

rJfi 

( 7 ) 


ds „ ^, 

— = -VsH X s. 
dt 


This allowed Breiter et ah (2005) to construct an efficient Lie-Poisson integrator for a fast 


propagation of the secular evolution of planetary spins. In Section we will use the leap¬ 
frog variant of the Hamiltonian’s LP2 splitting from Breiter et ah | (|2005). To propagate 


s through a single integration step, Breiter et al. (2005) method requires that the orbital 


semimajor axis, eccentricity and c are provided at times corresponding to the beginning and 
end of the step. These values can be supplied from an analytic model of orbit evolution, or 
be directly obtained from a previous numerical integrations of orbits where a, e and c were 
recorded with a conveniently dense sampling. 


2.2. Parametrization using obliquity and precession angle 


The Hamiltonian formulation in Eq. ([^ allows us to introduce several important con¬ 
cepts of the Cassini dynamics. A long tradition in astronomy is to represent s with obliquity 
e and precession angle V’ such that s^ = (sinesin-0, sine:cos"0, cose). The benefit of this 
parametrization is that the unit spin vector is expressed using only two variables. The 
drawback is that resulting equations are singular when e = 0. 


The conjugated momentum to '0 is = cose. The Hamiltonian is then (e.g., Laskar & 


Robutel 1993) 




-X'^ - 2CX 
2 


( 8 ) 















+\/l — X‘^ {A sin '0 + i3 cos -0). 


For the low-inclination orbits, C is negligible, while A and B are expanded in the Fourier 
series with the same frequency terms as those appearing in itself. A model of fundamental 
importance, introduced by G. Colombo (Colomb^|1966; Henrard fc Murigande ||1987 Ward 


& Hamilton 2004), is obtained when only one Fourier term in A and B is considered. 


This Colombo model obviously serves only as an approximation of the complete spin axis 
evolution, since all other Fourier terms in A and B act as a perturbation. Nevertheless, the 
Colombo model allows us to introduce several important concepts that are the basis of the 
discussion in Section [H 


In the Colombo model, the orbital inclination is hxed and the node precesses with a 
constant frequency. Put in a compact way, C = Aexp[*(st -f 0)] is the single Fourier term, 
and A = sin J/2 and VL = st + (p. Transformation to new canonical variables X' = —X 
and if = —{'ip + Q), and scaling with the nodal frequency s, results in a time-independent 
Hamiltonian 


nixpcp) = kX'^-cosIX' 

+ sin /Vl — cos ip, 


(9) 


where k = a/{2s) is a dimensionless parameter. Note that the orbit-plane angle (p is mea¬ 
sured from a reference direction that is 90° ahead of the ascending node. The general 
structure of the phase flow of solutions = C, with C constant, derives from the 

location of the stationary points. Depending on the parameter values (k, /), there are either 
two (|k| < K*) or four (|fi:| > k*) such stationary solutions (called the Cassini states). The 
critical value of n reads (e.g., Henrard & Murigande 1987|) 


K*(/) = - (sin2/3 j + 


cos 


2/3 


I) 


3/2 


( 10 ) 


Therefore, for small /, k* ~ The stationary solutions are located at (p = 0° or (p = 180° 
meridians in the orbital frame, and have obliquity values given by a transcendental equation 


K sin 2e = — sin (e =|= /), 


( 11 ) 


with the upper sign for ip = 0°, and the lower sign for ip = 180°. 

In the present work, we are mainly interested in the Cassini state C 2 located at (p = 0°. 
For Jupiter, the C 2 state related to frequency s = sy is subcritical since |k| < 1 for all 
estimates of the Jupiter’s precession constant found in the literature. Only two Cassini states 
exist in this regime, and s must circulate about C 2 . In the case of Saturn, |fi:| > k* ~ 1 
for s = Sg, four Cassini states exist in this situation, and s was suggested to librate in the 
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resonant zone abont C 2 . The confignration of vectors in C 2 can be inferred from Eq. 0. If 
the inclination is signihcantly smaller than the obliquity e, we have that a cose ~ —s. Since 
the term on the left hand side of this relation is the precession frequency of the planet’s spin 
Eq. [^, we hnd that the spin and orbit vectors will co-precess with the same rate about 


see 


the normal vector to the inertial frame. Small resonant librations about C 2 would reveal 
themselves by small departures of the spin vector from this ideal state. The maximum width 
Ae of the resonant zone in obliquity at the (p = 0 ° meridian can be obtained using an analytic 


formula (e.g., Henrard & Murigande 1987) 


Ae 


sm — = 


1 /sin 2 / 
IkI V sin 2 ^ 4 ’ 


( 12 ) 


where £4 is the obliquity of the unstable Cassini state C 4 (a solution of Eq. (11) at (p = 180° 
having the intermediate value of the obliquity). 

Figure [1} top panel, shows how the location of the Cassini states and the resonance 
width Ae depend on k, which is the fundamental parameter that changes during planetary 
migration. For sake of this example we assumed orbital inclination I = 0.5° (note that the 
overall structure remains similar for even smaller inclination values considered in the next 
section, but would be only less apparent in the Figure). The Ci and C4 stationary solution 


bifurcate when k = k* at a non-zero critical obliquity value e* = atan(tan^/^ I) (e.g., Henrard 


& Murigande 1987 Ward & Hamilton 2004). Note that As is signihcant in spite of a very 


small value of the inclination, which manifests through its dependence on a square root of 
sin 2 / in ( 12 ). The bottom panels show examples of the phase portraits 77 (X', 9 ?) = C for 
both sub-critical |fi:| < k* and super-critical |k| > k* cases. 


3. Results 

We now turn our attention to the evolution of Jupiter’s and Saturn’s obliquities dur¬ 
ing planetary migration. We hrst discuss the orbital evolution of planets in the instabil¬ 
ity/migration simulations of NM12 (Section 3.1). We then parametrize the planetary mi¬ 
gration before (stage 1 ) and after the instability (stage 2 ), and use it to study the effects on 
Jupiter’s and Saturn’s obliquities. The two stages are considered separately in Sections 3.2 
and 3.3. 
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3.1. The orbital evolution of giant planets 


NM12 reported the results of nearly 10^ numerical integrations of planetary instability, 
starting from hundreds of different initial configurations of planets that were obtained from 
previous hydrodynamical and A^-body calculations. The initial configurations with the 3:2 
Jupiter-Saturn mean motion resonance were given special attention, because Jupiter and 
Saturn, radially migrating in the gas disk before its dispersal, should have become trapped 


into their mutual 3:2 resonance (e.g., Masset & Snellgrove 2001 Morbidelli & Crida 2007 


Pierens & Nelson 2008). They considered cases with four, five and six initial planets. 


where the additional planets were placed onto resonant orbits between Saturn and the inner 
ice giant, or beyond the orbit of the outer ice giant. The integrations included the effects 
of the transplanetary planetesimal disk. NM12 experimented with different disk masses, 
density profiles, instability triggers, etc., in an attempt to find solutions that satisfy several 
constraints, such as the orbital configuration of the outer planets, survival of the terrestrial 
planet, and the distribution of orbits in the asteroid belt. 

NM12 found the dynamical evolution in the four planet case was typically too violent 
if Jupiter and Saturn start in the 3:2 resonance, leading to ejection of at least one ice giant 
from the Solar System. Planet ejection can be avoided if the mass of the transplanetary disk 
of planetesimals was large (Mdisk ^ 50 MEarth, where MEarth is the Earth mass), but such a 
massive disk would lead to excessive dynamical damping (e.g., the outer planet orbits become 
more circular then they are in reality) and to smooth migration that violates constraints from 
the survival of the terrestrial planets, and the asteroid belt. Better results were obtained 
when the Solar System was assumed to have five giant planets initially and one ice giant, with 
mass comparable to that of Uranus or Neptune, was ejected into interstellar space by Jupiter 
(Nesvorny 2011| Batygin et al. 2012). The best results were obtained when the ejected 
planet was placed into the external 3:2 or 4:3 resonances with Saturn and Mdisk — 20 MEarth- 
The range of possible outcomes was rather broad in this case, indicating that the present 
Solar System is neither a typical nor expected result for a given initial state. 


The most relevant feature of the NM12 models for this work is that the planetary 
migration happens in two stages (see Fig. [^. During the first stage, that is before the 
instability happens, Neptune migrates into the outer disk at ~ 20 — 30 an. The migration 
is relatively fast during this stage, because the outer disk still has a relatively large mass. 
We analyzed several simulations from NM12 and found that Neptune’s migration can be 
approximately described by an exponential with the e-folding timescale r ~ 10 Myr for 
Mdisk = 20 MEarth and r ~ 20 Myr for Mdisk = 15 MEarth- The instability typically happens in 
the NM12 models when Neptune reaches ~ 28 an. The main characteristic of the instability 
is that planetary encounters occur, mainly between the extra ice giant and all other planets. 
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The instability typically lasts ~ 10^ years and terminates when the extra ice giant is ejected 
from the solar system by Jupiter. The second stage of migration starts after that. The 
migration of Neptune is much slower during this period, because the outer disk is now very 
much depleted. From simulations in NM12 we find that r ~ 30 Myr for Mdisk = 20 MEarth 
and r ~ 50 Myr for Mdisk = 15 MEarth- Moreover, rather then being precisely exponential, 
the migration slows down relative to an exponential with fixed r, such as, effectively, the 
very late stages have larger r values (r ~ 100 Myr) than the ones immediately following 
the instability. Uranus accompanies the migration of Neptune on timescales similar to those 
mentioned above. 


The frequencies S 7 and Sg, which are the most relevant for this work, are initially some¬ 
what higher than the precession constants of Jupiter and Saturn, mainly because of the 
torques from the outer disk. The extra term from the third ice giant initially located at 
~ 10 an has much faster frequency then the precession constants and does not interfere 
with the obliquity of Jupiter and Saturn during the subsequent evolution. The sy and sg 
frequencies slowly decrease during both stages. Their e-folding timescales may slightly differ 
from the migration e-folding timescales mentioned above, due to the non-linearity of the 
dependence of the secular frequencies on the semimajor axis of planets. Our tests show that 
they are about (90 — 95)% of the e-folding timescales of planetary semimajor axes. 


From analyzing the behavior of frequencies in different simulations we found that Sg 
should cross the value of Jupiter’s precession constant during the hrst migration stage, 
that is before the instability. The main characteristic of this crossing is that the planetary 
orbits are very nearly coplanar during this stage. The amplitude / in Eq. 0 should thus 
be very small. It is not known exactly, however, how small. In Section 3.2, we consider 
amplitudes down to 0.005° (about 10 times smaller than the present value of /58, where /58 
is the amplitude of the 8 th frequency term in Jupiter’s orbit), and show that the effects on 
Jupiter’s obliquity are negligible if the amplitudes were even lower. The second characteristic 
of the hrst migration stage is that the evolution of sg happens on a characteristic timescale 
of ~ 10 — 20 Myr. Since the total change of Sg during this interval is several arcsec yr“^, and 
the hrst stage typically lasts ~ 10 — 20 Myr, the average rate of change is very roughly, as 
an order of magnitude estimate, dsg/dt ~ 0.1 arcsec yr“^ Myr“^. The actual value of dsg/dt 
during crossing depends on several unknowns, including when exactly the crossing happens 
during the hrst stage. Also, the changes of sg could have been slower if the hrst stage lasted 
longer than in the NM12 simulations, as required if the instability occurred at the time of 
the Late Heavy Bombardment 


fe-g. 


Gomes et ah 


2005). In Section 3.2, we will consider 
values in the range 0.005 < dsg/dt < 0.05 arcsec yr“^ Myr“^, and show that the obliquity 
of Jupiter cannot be pumped up to its current value if dsg/df > 0.05 arcsec yr“^ Myr“^ 
(assuming that Jsg < 0.05°). 
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Interesting effects on obliqnities shonld happen during the second migration stage. First, 
the Sg frequency reaches the value of the precession constant of Saturn. There are several 
differences with respect to the sg crossing of Jupiter’s precession constant during the hrst 
stage (discussed above). The orbital inclinations of planets were presumably excited to their 
current values during the instability. Therefore, the amplitude /eg should be comparable to 
its current value, Jgg ~ 0.064°, during the second stage. We see this happening in the NM12 
simulations. First, there is a brief period during the instability, when the inclinations of all 
planets are excited by encounters with the ejected ice giant. The inclination of Neptune is 
modest, at most ~ 2 °, and is rather quickly damped by the planetesimal disk. Also, the 
invariant plane of the solar system changes by ~ 0.5° when the third ice giant is ejected 
during the instability. The hnal inclinations are of this order. The current amplitudes are 
Jsg ~ 0.066° (sg term in Jupiter’s orbit) and Jgg — 0.064° (sg term in Saturn’s orbit) (see 
e.g. ^ " 


Laskar 1988). 


Another difference with respect the hrst stage is that the evolution of sg is much slower 
during the second stage. If, as indicated by the NM12 integration, sg changes by ~ 1 arc- 
sec yr“^ in 100 Myr, then the average rate of change is very roughly dsg/dt ~ 0.01 arc- 
sec yr“^ Myr“^. The actual rate of change can be considerably lower than this during the 
very late times, when the effective r was lower than during the initial stages. Finally, during 
the very last gasp of migration, the S 7 frequency should have approached the precession con¬ 
stant of Jupiter. We study this case in an adiabatic approximation when the rate of change 
of S 7 is much slower than any other relevant timescale. We hud that the present obliquity of 
Jupiter can be excited by the interaction with the S 7 term only if the precession constant of 


Jupiter is somewhat larger than inferred by Helled et ah (2011), in accord with the results 


of Ward & Canup (2006). 


3.2. The effects on Jupiter’s obliquity during stage 1 

Since sg remains larger than as during the hrst stage, we do not expect any important 
effects on Saturn’s obliquity during this stage. If Saturn’s obliquity was low initially, it 
should have remained low in all times before the instability. We therefore focus on the case 
of Jupiter in this section. From the analysis of the NM12 numerical simulation in Section 3.1, 
we infer that the Sg frequency crossed oj during the hrst stage. The values of dsg/dt and /sg 
during crossing are not known exactly from the NM12 simulations, because they depend on 
details of the initial conditions. We therefore consider a range of values and determine how 
Jupiter’s obliquity excitation depends on them. The results can be used to constrain future 
simulations of the planetary instability/migration. 
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We consider Colombo’s model with only one Fonrier term in ( of Jupiter, namely that of 
the Sg frequency]^ The inclination term in Jupiter’s orbit was treated as a free parameter. 
The range of values was set to be between zero and roughly the current value of ~ 0.066°. As 
we discussed in Section 3.1, it is reasonable to assume that I^s was smaller than the current 
value, because the orbital inclination of planets should have been low in times before the 
instability. The value of aj was obtained by rescaling the present value to the the semimajor 
axis of Jupiter before the instability (oj ~ 5.45 an). To do so we used Eq. Q and assumed 
aj oc aj"^. No additional modeling of possible past changes of ctj, for instance due to satellite 
system evolution or planetary contraction, was implemented. The sg frequency was slowly 
decreased from a value larger than aj to a value smaller than aj. 

Motivated by the numerical simulations of the instability discussed in Section 3.1, we 
assumed the initial sg value of —4 arcsec yr“^ and let it decrease to —1.2 arcsec yr“^ by 
the end of each test. The rate of change, dsg/dt, was treated as a free parameter. The 
integrations were carried for several tens of Myr for the highest assumed rates and up to 
several hundreds of Myr for the lowest rates. We recorded Jupiter’s obliquity during the 
last 5 Myr of each run, and computed the mean value egn- Jupiter’s initial spin axis was 
oriented toward the pole of the Laplacian plane. [The value of £fin reported in Fig. |^was 
averaged over all possible phases of the initial spin axis on the ip) = C level curve, 

with C dehned by s oriented toward the pole of the Laplacian plane] 

Figure shows the results. For most parameter combinations shown here the Sg res¬ 
onance swept over aj without having the ability to capture Jupiter’s spin vector in the 
resonance. This happened because dsg/df was relatively large and I^g was relatively small, 
thus implying that the Sg frequency crossed the resonant zone in a time interval that was 
shorter than the libration period. Captures occurred only in extreme cases (largest I^g and 
smallest dsg/dt). These cases ended up generating very large obliquity values of Jupiter and 
are clearly implausible. The plausible values of dsg/df and /58 correspond to the cases where 
Jupiter’s obliquity was not excited at all, thus assuming that Jupiter obtained its present 
obliquity later, or was excited by up to ~ 3°. To obtain ej ~ 3°, dsg/dt and I^g would need 
to have values along the bold line labeled 3 in Fig. which extends diagonally in dsg/df 
and Jsg space. An example of a case, where the obliquity of Jupiter was excited to near ~ 3° 
value, is shown in Fig. 


^We found that adding higher frequency terms, such as sq and/or sj, into our simulation does not change 
results. 
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3.3. Behavior of obliquities during stage 2 


We now turn our attention to the second stage, when the migration slowed down and 
the obliquities of Jupiter and Saturn should have suffered additional perturbations. At the 
beginning of stage, that is just after the time of the instability, the Sg frequency is already 
lower than oj, but still higher than as, while the sj frequency is higher than oj. Since the 
sj and ss frequencies are slowly decreasing during the second stage, a possibility arises that 
Jupiter’s obliquity was (slightly) excited when sj approached aj (e.g.. Ward fc Canup [2006 ), 
and that Saturn’s obliquity was strongly excited by capture into the spin-orbit resonance 
with sg (e.g.. 


Ward & Hamilton 2004; Hamilton & Ward 2004; Bone et al. 2009). 


3.3.1. Jupiter 


Ward & Canup 


(2006) suggested a possibility that Jupiter’s present obliquity may be 
explained by the proximity of aj to the current value of the Sy frequency. They showed 
that, if K = aj/( 2 s 7 ) is sufficiently close to the critical value from Eq. (10), namely ~ ^ 
for small inclinations, the obliquity of the Cassini state C 2 may be significant. Thus, as 
sy adiabatically approached to aj, Jupiter’s obliquity may have been excited along. As a 
supportive argument for this scenario they pointed out that 93 J ~ 0°, where the Cassini state 
2 is located. 


To test this possibility we run a suite of simulations, assuming an exponential conver¬ 
gence to the current value sy ~ —2.985 arcsec yr“h Specifically, we set sy(t) = sy(0) -|- [sy — 
sy(0)] exp(—t/r), where r and sy(0) are parameters. The initial value sy(0) at the beginning 
of the second stage was obtained from the numerical simulation discussed in Section 3.1. Here 
we chose to use sy(0) = —3.5 arcsec yr“^, however we verified that the results are insensitive 
to this choice. The e-folding timescale r depends on how slow or fast planets migrate. Given 
that the planetary migration is slow during the second stage, we chose r = 100 — 200 Myr. 
This assures that the approach of sy(t) to aj is adiabatical. The amplitude Jsy is assumed 
to be constant and equal to its current value (Jsy ~ 0.055°). 

Two additional parameters need to be specified: (i) Jupiter’s initial spin state, and (ii) 
aj. As for (i), the results discussed in Section 3.2 indicate that Jupiter’s obliquity may 
have remained near zero during the first stage, if Jsg was too small and/or dsg/dt was too 
fast, or could have been potentially excited to ~ 3°, if Jsg and dsg/df combined in the right 
way. Therefore, here we treat the obliquity of Jupiter at the beginning of stage 2 as a 
free parameter. As for (ii), as we discussed in Section 1, the present value of aj somewhat 
uncertain. We therefore performed various simulations, where aj takes on a number of 
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different values between 2.75 arcsec yr ^ and 2.95 arcsec yr A similar approach has also 
been adopted by Ward & Canup ( 2006| . 


Figure [^reports the results. The top panel shows how the obliquity £c 2 of the Cassini 
state C 2 depends on the assumed value of aj. This is calculated from Eq. ( 0 . The trend is 
that ec 2 increases with aj, because the larger values of aj correspond to a situation where 
the system is closer to the exact resonance with sj. If aj < 2.8 arcsec yr“^, as inferred from 


models in Helled et al. (2011), ec 2 is too small to signihcantly contribute to ej. This case 


would imply that Jupiter’s present obliquity had to be acquired during the earlier stages 
and is possibly related to the non-adiabatic sg resonance crossing discussed in Section 3.2. 
If, on the other hand, aj ~ 2.92 — 2.94 arcsec yr“^, Jupiter would owe its present obliquity 
to the proximity between aj and S 7 . This would imply that the obliquity excitation during 
the hrst stage of planetary migration must have been minimal. Figures and [^express the 
joint constraint on the planetary migration also for the intermediate cases, where the present 
obliquity of Jupiter arose as a combination of both effects discussed here. 


Figure illustrates the two limiting cases discussed above. In panel (a), we assumed 
that the parameters during the hrst migration stage were such that the obliquity of Jupiter 
was excited to its current value during the sg crossing (such as shown in Fig. |^. Also, we set 
aj = 2.77 arcsec yr“^, corresponding to the best theoretical value of Helled et al. (2011). 
The Cassini state C 2 corresponding to the s^ term is only slightly displaced from the center 
of the plot, and does not signihcantly contribute to the present obliquity value. In panel (b) 
of Fig. we set aj = 2.93 arcsec yr“F This value implies that ec 2 — 2.6°. The present 
obliquity of Jupiter would then be in large part due to the “forced” term arising from the 
proximity of S 7 . The initial excitation of Jupiter’s obliquity during the hrst stage would have 
to be minor in this case. 


Helled et al. (2011 


3.3.2. Saturn 


In the case of Saturn, all action is expected to take place during the second stage of plan¬ 
etary migration. Insights gleaned from the numerical simulations discussed in Section 3.1 
show that the Sg frequency should have very slowly approached as, thus providing a concep¬ 


tual basis for capture of Saturn’s spin vector in a resonance with sg (e.g.. Ward & Hamilton 


2004; Hamilton & Ward 2004 Bone et al. 2009). To study this possibility, we assume 


that the leg amplitude was excited to its current value during the instability, and remained 
nearly constant during the second stage of planetary migration. This choice is motivated by 


the NM12 simulations, where the inclination of Neptune is never too large. Note that Bone 


et al. (2009) investigated the opposite case where Neptune’s inclination was substantially 
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excited during the instability and remained high when the resonance with sg occurred. This 
type of strong inclination excitation does not happen in the NM12 models. 


Here we assume that the planetary migration was very slow during the second stage 
and parametrize sgit) as sgit) = S8(0) + [sg — S8(0)] exp(—t/r) with r > 80 Myr. The initial 
frequency value at the beginning of stage 2, Sg(0), is estimated from the NM12 simulations. 
We find that S8(0) ~ —1.3 arcsec yr“^, and use this value to set up the evolution of Sgit). 
We also assume a range of as values. This has the following signihcance. As already pointed 


out by Boue et al. (2009), the best-modeling values of as from Helled et ah (2009) are 


not compatible with a resonant location of Saturn’s spin axis. This is because the Cassini 
state C 2 would be moved to a signihcantly larger obliquity value (> 34°). So these values 
of as would imply that Saturn’s spin circulates about the Cassini state Ci. On the other 
hand, the significant obliquity of Saturn requires an increase when the sg value was crossing 
as value, as schematically shown in the left panel of Fig. for Jupiter’s obliquity during the 


phase 1. Boue et al. (2009) tested this scenario using numerical simulations and found it 


extremely unlikely: initial data of an insignificant measure have led this way to the current 
spin state of Saturn. Indeed, here we recover the same result with a less extensive set of 
numerical simulations. 


Given the arguments discussed above we therefore tend to believe that the precession 


constant of Saturn may be somewhat smaller than the one determined by Helled et al. 


(2009). For instance, R. A. Jacobson (personal communication) determined the Saturn 
precession from the Saturn’s ring observations. The mean precession rate obtained by him 
is 0.725 arcsec yr“^ (formal uncertainty of about 6%). This value would indicate as in the 
range between 0.769 arcsec yr“^ and 0.864 arcsec yr“^. Both Ward & Hamilton (|2004) and 


Boue et al. (2009) report other observational constraints of Saturn’s pole precession that 


have comparably large uncertainty. We therefore sampled a larger interval of the as values 
to make sure that all interesting possibilities are accounted for. 


Our numerical simulations thus spanned a grid of two parameters: (i) as discussed 
above, and (ii) r, the e-folding timescale of the Sg frequency that slowly changes due to 
residual migration of Neptune and depletion of the outer disk. The amplitude related to the 
Sg term in the inclination vector ( of Saturn is kept constant, namely I^g = 0.064°. To keep 
number of tested free parameters low, we assumed initial orientation of Saturn’s spin axis s 
to be near the pole of the invariable plane. Specifically, we set its obliquity to 0.1° in the 
reference frame of the Sg-frequency Fourier term in (. To prevent fluke results, we sampled 
36 values of longitude (p of the Saturn’s pole in the same reference frame, incrementing it 
from 0° by 10°. Each of the simulations covered a 1 Gyr timespan. We recorded Saturn’s 
pole orientation during the last 150 Myr time interval. We specihcally analysed if it passes 
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close to the current location of Saturn’s pole, namely £§ — 27.4° and <p8 ~ —31.4° in the 
Sg-frequency reference frame (see Table 2 in Ward & Hamilton 2004). A numerical run 
was considered successful, if the simulated Saturn’s pole passed through a box of ±0.2° in 
obliquity e and ±3° in longitude (p around the planet’s values (eg, ips) during the recorded 
150 Myr time interval. Note that the libration period of Saturn’s pole around Cassini state 
C 2 , if captured in the spin-orbit resonance, is ~ (50 — 100) Myr, depending on the libration 
amplitude. This set our requirement for the timespan over which we monitored Saturn’s 
pole position. 


Figure shows the results from this suite of runs. The shaded region shows correlated 
as-T pairs that provided successful match to the Saturn’s pole position. We note that no 
successful solutions were obtained for as > 0.812 arcsec yr“^ and all successful solutions 
correspond to the capture in the resonance zone around the Cassini state C 2 . The absence of 
low-probability solutions in which Saturn’s pole would circulate about the Cassini state Ci 
for larger as may be related to the limited number of simulations performed. No solutions 
were also obtained for r > 215 Myr. This is because for such long e-folding timescales the 
resonant capture process would be strictly adiabatic and the simulated spin would not meet 


the condition of at least 30° libration amplitude (see discussion in Ward & Hamilton 2004 


Hamilton fc Ward~||2004 ). The area occupied by successful solutions splits into two branches 
for as — 0.78 arcsec yr“^. This is because the obliquity of the Cassini state C 2 is ec 2 — 27.4° 
for the critical value of as, and the solutions have the minimum libration amplitude of ~ 31°. 


Figure shows two evolution paths of Saturn’s spin vector obtained in two different 
simulations. As mentioned above, in both cases Saturn’s spin state was captured into the 
spin-orbit resonance Sg, and remained in the resonance for the full length of the simulation. 
The hnal states of both simulations are a good proxy for the Saturn’s present spin state. 
These two cases differ from each other principally because the path in panel (b) shows li- 
brations with a larger amplitude than the path in panel (a). Note some of the solutions, 
such as (b) here, may attain a signihcant librations amplitude. This is because with the 
corresponding values of r ~ 100 Myr the evolution is not adiabatic and the librations am¬ 
plitude is excited immediately after capture. Therefore, the complicated evolution histories 
proposed in Hamilton & Ward (2004) may not be needed. 


4. Conclusions 

We studied the behavior of Jupiter’s and Saturn’s obliquities in models of planetary 
instability and migration that were informed from NM12. Rather then investigating a few 
specihc cases directly from NM12, we considered the general concept of a two-stage migration 









18 


from NM12, and studied a broad range of relevant parameter values. We found that, in 
general, the two stage migration provides the right framework for an adequate excitation 
of Jupiter’s and Saturn’s obliquities. Moreover, we found that certain conditions must be 
satisfied during the first and second stages of migration, if the final obliquity values are to 
match the present obliquities of these planets. 

Our results indicate that Jupiter spin axis could have been tilted either when (i) the 
ss frequency swept over aj during the hrst migration stage (that is before the instability 
happened), or when (ii) the sy frequency approached aj at the end of planetary migration. 
For (i) to work, the crossing of sg must be fast, such that the capture into the resonance 
does not happen, but not too fast, such that some excitation is generated by the resonance 
crossing. To obtain full ej = 3.1° during this stage, the rate of change of Sg during crossing, 
dss/df, must be smaller than 0.05 arcsec yr“^ Myr“^ (assuming that I^g < 0.05°). Since 
the evolution of sg mainly relates to the radial migration of Neptune and dispersal of the 
outer disk of planetesimals, this result implies that both these processes would need to occur 
relatively slowly. More specifically, parameters dsg/dt and I^g would have to have values 
along a diagonal line in the (dsg/df, I^g) plane with larger values of I^g requiring larger 
values of dss/df (see Fig. [^. Any model of planetary instability/migration can be tested 
against this constraint. The models where dss/dt is too slow and/or I^g is too large, as 
specihed in Fig. can be ruled out, because Jupiter’s obliquity would be excited too much 
by the Sg crossing. 


Not much excitation of ej is expected during the sg crossing if dss/dt was relatively fast 
and/or if I^g was only a very small fraction of its current value. If that is the case, Jupiter’s 
obliquity would probably need to be excited during the very last stages of migration by 
(ii). For that to work, Jupiter’s precession constant aj would need to be ~ 2.95 arcsec 
yr“^ (assuming £j = 0 initially), which is a value that is signihcantly larger than the one 


estimated by Helled et ah (2011). This means that Helled’s model would need to be adjusted 


to £t within this picture. It is also possible, however, that Jupiter’s present obliquity was 
contributed partly by (i) and partly by (ii). If so. Figs. and express the joint constraint 
on dsg/dt and I^g during the hrst stage, and aj. 

As for Saturn, our results indicate that the capture into the spin-orbit resonance with 


Sg (Ward & Hamilton 2004; Hamilton & Ward 2004) is indeed possible during the late 


stages of planetary migration, assuming that the migration rate was slow enough. The exact 
constraint on the slowness of migration depends on Igg, which in turn depends how much 


Neptune’s inclination was excited by the instability and how long it remained elevated (Bone 


et ah 2009). Since in the NM12 models, Neptune’s orbital inclination is never large, we 


have good reasons to believe that Igg was comparable to its current value when the crossing 
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of Sg occurred. Thus, using Igg — 0.064°, we find that the e-folding migration timescale r 
would need to be r > 100 Myr. If r > 200 Myr, however, the capture in the Sg resonance 
would be strictly adiabatic. This would imply, if £s was negligible before capture, that the 
resonant state should have a very small libration amplitude (see Ward & Hamilton 2004 


Hamilton & Ward 2004). It would then be difficult to explain the current orientation of 
Saturn’s spin axis, which indicates that the libration amplitude should be at least 30°. A 
more satisfactory solution, however, can be obtained for r ~ 100 — 150 Myr, in which case 
the capture into the resonance was not strictly adiabatic. In this case the > 30° libration 
amplitude is obtained during capture. 

While the capture conditions pose a strong constraint on the timescale of Neptune’s ra¬ 
dial migration, as discussed above, and additional constraint on Saturn’s precession constant 
as derives from the present obliquity of Saturn. This is because, again assuming that Saturn 
spin vector is in the resonance with sg today, the present obliquity of Saturn implies that the 
Cassini state C 2 of this resonance would have to be located at £ ~ (28 — 30)°. This would 
require that as — 0.78 — 0.80 arcsec yr“^. The value derived by Helled et ah (2009) is larger, 
~ 0.8445 arcsec yr“^, and clearly incompatible with this assumption. Direct measurements 
of the mean precession rate of Saturn’s spin axis suggest that as — 0.81 arcsec yr“^, which is 
still slightly larger than the range given above, but the uncertainty interval of this estimate 
includes values below 0.8 arcsec yr“^ (R. A. Jacobson, personal communication). Figuring 
out the exact value of Saturn’s precession constant will therefore be important. Once as is 
known. Fig. could be used to precisely constrain the timescale of planetary migration. 


We are grateful to Robert Jacobson for sharing his latest estimate of the precession 
rate of Saturn’s spin axis. We also thank Tristan Guillot and Alessandro Morbidelli for 
discussions and pointing to us the future capability of Juno mission to constrain Jupiter’s 
precession constant. The work of DV was supported by the Czech Science Foundation (grant 
GA13-01308S). The work of DN was supported by NASA’s Origin of Solar Systems (OSS) 
program. 
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Fig. 1.— Top panel: Parametric dependence of Cassini state Ci, C 2 and C 4 obliquity 
(ordinate) on the frequency-ratio parameter k (abscissa) in the Colombo model (the C 3 
equilibirium has obliquity larger than 90° and it is not relevant for our discussion). The 
orbital inclination I has been set to 0.5° value. The dashed line indicates the critical value 
—K* at which Ci and C 4 bifurcate (Eq. 10). The gray area for \k\ > shows maximum width 
of the resonant zone around C 2 . Bottom panels: Examples of phase portraits 'H(X', (p) = C 
for two values of k: (a) k, = —0.4 on the left, and (b) k, = —0.6 on the right. We use 
coordinates x = sine cos 99 and y = sine sin 99 with the origin at the north pole e = 0. The 
gray symbols show location of the Cassini equilibria and the curves are isolines 'H{X', p) = C 
for suitably chosen C values. The bold line in (b) is the separatrix of the resonant zone around 
C 2 stationary solution. 
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Fig. 2.— An example of planetary migration and instability from Nesvorny fc Morbidelli 


(2012). The plot shows the evolntion of the semimajor axis (bold line), and the perihelion 
and aphelion distances (thin lines) of the giant planets. The initial orbits of Jnpiter, Satnrn 
and the inner ice giant were placed in the 3:2 resonant chain. The semimajor axes of the two 
onter ice giants were set to be 16 an and 22 an. The trans-Neptnnian disk of planetesimals 
(not shown here) was resolved by 10,000 eqnal-mass particles. The disk originally extended 
from 23.5 an to 29 an and had the total mass of 20 MEarth- The instability happened at 
t ~ 5.6 Myr after the start of the simnlation. The third ice giant was snbseqnently ejected 
from the Solar system. Note that the migration rate before the instability (phase 1) is 
signihcantly larger than the migration rate after the instability (phase 2). 
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Fig. 3.— The final obliquity egn of Jupiter obtained in our integrations of crossing of the 
Ss spin-orbit resonance. Jupiter’s obliquity is shown as a function of two parameters: (i) 
the amplitude I^s of the Fourier term in Jupiter’s ( variable (see Section 2), and (ii) the 
rate dsg/dt at the time when sg became equal to aj. The gray shading indicates the final 
obliquity (the scale bar on the right shows the corresponding numerical value in degrees). 
The three bold isolines correspond to ej = 1°, 2° and 3° (see the labels). 



25 



Fig. 4.— An example demonstrating the effect of ss frequency sweeping over aj. Here we set 
Jss = 0.02° and dsg/dt = 0.01 arcsec yr“^ Myr“F The left panel shows Jupiter’s obliquity as 
a function of time. Obliquity £j increases to ~ 2.7° during the resonance crossing. The width 
of the Cassini resonance is small because I^s is small, and the assumed rate dss/dt is too 
large in this case to lead to capture. The right panel shows the spin axis evolution projected 
onto the {x,y) plane, where x = sine cos 99 and y = sine sin 99 . The Cassini state C 2 drifts 
along the x-axis during the integration (as indicated by the gray arrow), and reaches very 
large obliquity values. 
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Fig. 5.— Final obliquity of Jupiter resulting from an adiabatic approach of the Sy-frequency 
towards aj (the obliquity has been computed in the reference frame associated with this 
frequency term in Q. While today’s value of sy is known fairly accurately, aj has a signihcant 
uncertainty. We therefore treat aj as a free parameter. The initial obliquity of Jupiter is 
also treated as a free parameter, because it depends on the effects during the hrst migration 
stage (see Fig. [^. The key to the shading scale is provided by the vertical bar on the right 
(white region corresponds to the hnal maximum obliquity smaller than 2.5°, incompatible 
with Jupiter’s value). The bold curve corresponds to the 3.45° isoline, estimated obliquity 
value today in this frame (e.g.. Ward & Canup 2006|. The arrows indicate parametric 


location of the two examples shown on Fig. The upper panel shows the obliquity value of 
the Cassini state C 2 as a function of aj. 
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Fig. 6.— Two examples of Jupiter’s spin state evolution during the 1 Gyr-long time interval 
after the planetary instability. The planet pole is shown in the Cartesian coordinates x = 
sine cos (p and y = sine sin (p. The arrow shows evolution of the Cassini state C 2 over the 
interval of time covered by the simulation. The gray star shows the current location of 
Jupiter’s pole in these coordinates. Cray dots show the pole position output every 5 kyr 
during the simulation, the black circles highlight the hrst and the last 30 Myr of the evolution. 
Two different parametric combination along the bold curve in Fig. [^were chosen: (i) aj = 
2.77 arcsec yr“^ and £ini = 3.1° in the left panel (a), and (ii) aj = 2.93 arcsec yr“^ and 
£ini = 1.3° in the right panel (b). 
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Fig. 7.— Bottom panel: Distribution of solutions successfully matching Saturn’s spin state 
in the parametric plane as (abscissa) vs r (ordinate), the e-folding time of sg-frequency 
evolution during the phase 2. No successful solutions were obtained in the white region of 
the plot. The darker the gray-scale in the given bin, the more robust the solution is. The 
maximum value 36 corresponds to 36 sampled initial conditions of the simulations for each 
(as, r) pair (see the side bar). When as ~ 0.78 arcsec yr“^, the obliquity of the corresponding 
to Cassini state C 2 is ~ 27.4°. Therefore the capture solution corresponds to the minimum 
needed libration amplitude of Saturn’s spin in the Sg reference frame. From this value the 
larger libration-amplitude solutions bifurcate towards smaller/larger as values. The arrows 
indicate parametric location of two particular examples shown on panels (a) and (b) of Fig. 
Top panel: Obliquity of the Cassini state C 2 (solid line) and maximum width of the associated 
resonant zone (gray area) as a function of as- We assume inclination / = Jgg = 0.064° and 
terminal value of the orbital precession frequency Sg ~ —0.692 arcsec yr“^. 
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Fig. 8.— Two examples Saturn’s obliquity excitation by the capture and evolution in the 
resonance with sg. The spin axis s is projected onto the {x^y) plane, where x = sine:cos(p 
and y = sin £ sin (f. The reference frame used here is the one associated with the sg frequency 
term contributing to C, of Saturn. The gray arrow shows the evolution of the Cassini state C 2 
over the full length of the simulation (1 Gyr). In panel (a), we assumed that as = 0.785 arcsec 
yr“^, which means that the terminal obliquity of C 2 is at ~ 28.2°, and r = 150 Myr. In (b), 
we used as = 0.8 arcsec yr“^, which means that the terminal obliquity of C 2 is at ~ 30.3°, 
and r = 135 My. These combinations were also highlighted by arrows on Fig. The star in 
each panel shows the present orientation of the Saturn pole vector in the (x, y) coordinates 
(e.g.. Ward & Hamilton 2004, Table 2). 


















